"""
Allen-Cahn PDE implementation

The Allen-Cahn equation describes phase separation and interface dynamics:
    ∂c/∂t = M[γ∇²c + c - c³]

This is a gradient flow for the Cahn-Hilliard free energy with non-conserved dynamics.
"""

using ..Fields
using ..Grids

"""
    AllenCahnPDE <: AbstractPDE

Allen-Cahn equation for phase field modeling and interface dynamics.

The mathematical definition is:
    ∂c/∂t = M[γ∇²c + c - c³]

where c is the order parameter, M is the mobility, and γ controls the interface width.

# Parameters
- `interface_width`: Interface width parameter γ (default: 1.0)
- `mobility`: Mobility parameter M (default: 1.0)
- `bc`: Boundary conditions (default: "auto_periodic_neumann")

# Properties
- Second-order nonlinear PDE
- Energy decreasing: dE/dt ≤ 0
- Stable interfaces between c = ±1 phases
- Interface width ~ √γ
"""
struct AllenCahnPDE <: AbstractPDE
    interface_width::Float64  # γ parameter
    mobility::Float64         # M parameter  
    bc::String                # boundary conditions
    
    function AllenCahnPDE(; interface_width=1.0, mobility=1.0, bc="auto_periodic_neumann")
        new(interface_width, mobility, bc)
    end
end

function evolution_rate(pde::AllenCahnPDE, state::ScalarField, t::Float64)
    c = data(state)
    grid_obj = grid(state)
    
    # Diffusion term: γ∇²c
    c_laplace = laplacian(state)
    diffusion_term = pde.interface_width * c_laplace
    
    # Nonlinear potential term: c - c³ = c(1 - c²)
    potential_data = c - c.^3
    potential_term = ScalarField(potential_data, grid_obj)
    
    # Combined chemical potential: γ∇²c + c - c³
    chemical_potential = diffusion_term + potential_term
    
    # Apply mobility: ∂c/∂t = M * chemical_potential
    return pde.mobility * chemical_potential
end

"""
    free_energy(pde::AllenCahnPDE, state::ScalarField)

Compute the free energy functional:
    E = ∫[f(c) + (γ/2)|∇c|²] dV
where f(c) = (c² - 1)²/4
"""
function free_energy(pde::AllenCahnPDE, state::ScalarField)
    c = data(state)
    
    # Bulk free energy density: f(c) = (c² - 1)²/4
    bulk_energy = sum((c.^2 .- 1).^2) / 4
    
    # Gradient energy density: (γ/2)|∇c|²
    grad_c = gradient(state)
    grad_magnitude_sq = 0.0
    if isa(grad_c, VectorField)
        grad_data = data(grad_c)
        if ndims(grad_data) >= 2
            for i in 1:size(grad_data, 1)
                grad_magnitude_sq += sum(grad_data[i, :].^2)
            end
        else
            grad_magnitude_sq += sum(grad_data.^2)
        end
    end
    gradient_energy = pde.interface_width * grad_magnitude_sq / 2
    
    return bulk_energy + gradient_energy
end

"""
    get_interface_width(pde::AllenCahnPDE)

Get the characteristic interface width: w = √(2γ)
"""
get_interface_width(pde::AllenCahnPDE) = sqrt(2 * pde.interface_width)

"""
    is_energy_decreasing(pde::AllenCahnPDE)

Allen-Cahn dynamics decrease the free energy.
"""
is_energy_decreasing(pde::AllenCahnPDE) = true

"""
    expression(pde::AllenCahnPDE)

Return a string representation of the Allen-Cahn PDE.
"""
function expression(pde::AllenCahnPDE)
    terms = String[]
    
    # Diffusion term
    if pde.interface_width != 0
        γ_str = pde.interface_width == 1.0 ? "" : "$(pde.interface_width)"
        push!(terms, "$(γ_str)∇²c")
    end
    
    # Potential terms
    push!(terms, "c - c³")
    
    expr = join(terms, " + ")
    
    # Apply mobility
    if pde.mobility == 1.0
        return expr
    else
        return "$(pde.mobility)($(expr))"
    end
end

# Export the enhanced Allen-Cahn implementation
export AllenCahnPDE, free_energy, get_interface_width, is_energy_decreasing